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Abstract. We present results of a deep mid-infrared survey in the SSA13 field with the Infrared Space Observatory (ISO). 
In order to probe the near-infrared light at high redshifts, we surveyed the field with the broad band LW2 (5-8.5 yum) filter 
of the mid-infrared camera ISOCAM. Adopting a highly redundant imaging strategy for the 23 hour observation and carefully 
treating gradual changes in the detector responsivity caused by a very high rate of cosmic ray impacts, we succeeded in reaching 
an 80 % completeness limit of 16yuJy in the central 7 arcmin 2 region. Utilizing the signal-to-noise ratio map, we detected 65 
sources down to 6/jJy in the 16 arcmin 2 field. Integral galaxy number counts at 6.7 yum are then derived, reaching 1.3 x 10 4 
deg -2 at the faint limit with a slope of -1.6 between 13 pJy and 130/jJy. Integrating individual sources in this flux range, the 
resolved fraction of the extragalactic background light at 6.7 yum is estimated to be 0.56 nW irr 2 sr~' . These results, which reach 
a flux limit three times fainter than those in the Hubble Deep Fields, are in fairly good agreement with a model prediction by 
Franceschini et al. (1997). Finally, we discuss the relation of distant massive E/S0 galaxies to the faint 6.7 /an galaxy population. 

Key words. Galaxy: evolution - Infrared: galaxies - Cosmology: observations - Surveys - Galaxies: photometry - Methods: 
data analysis 



1. Introduction 

Explorations of the evolution of galaxies require representative 
galaxy samples at high redshifts. Intrinsically sensitive optical 
observations succeeded in picking up z ~ 3 sources utilizing 
characteristic breaks in galaxy spectra (Steidel et al. [1996 1. 
Another technique to search for high redshift galaxies is to 
observe at wavelengths long-ward of a global peak in spec- 
tral energy distributions (SEDs) of galaxies. A famous exam- 
ple is the recent discovery of distant dusty starbursting galax- 
ies in the Submillimeter and millimeter with the Submillimeter 
Common-User Bolometer Array SCUBA (Smail et al. 119971 
Barger et al. 119981 Hughes et al. 11998b and the Max-Planck 
Millimeter Bolometer array MAMBO (Bertoldi et al. 2000}. At 
submillimeter and millimeter wavelengths, contributions from 
foreground objects are small, since their dust emission has a 
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peak around 100 /im and decreases toward longer wavelengths. 
On the other hand, contributions from higher redshift objects 
become greater because of the K-correction brightening; the 
rest-frame 100/im dust peak shifts into the submillimeter or 
millimeter observing band with increasing redshift. As a result, 
the fraction of high redshift objects in the submillimeter or mil- 
limeter is greater than that at, say, 200 yum. 

The same technique can be applied to another striking peak 
around 1 yum originating from stellar photospheres. The ob- 
serving wavelength should be much longer than 1 //m to get 
the K-correction brightening. However, emission from PAH 
molecules and very hot dust should also be considered. Their 
contributions begin to increase at wavelengths long-ward of 
6 //m for local objects. We selected the broad band LW2 fil- 
ter at 6.7 //m in the mid-infrared camera ISOCAM (Cesarsky 
et al. I1996i as a good window between local stellar and dust 
emission. The fraction of high redshift objects would be larger 
at 6.7 //m than at 2.2 /mi. 

Previous field galaxy surveys at 6.7/im were carried out in 
the Hubble Deep Field (HDF) and the Lockman Hole (Serjeant 
et al. 119971 Taniguchi et al. 119971 . These found a few tens 
of 6.7 //m sources down to 30/dy (Goldschmidt et al. [1997 1, 
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some of which were confirmed by re-calibration of the HDF 
data (Desert et al. 119991 Aussel et al. I1999i . A more shallow 
(> 100 /(Jy) survey was conducted in the CFRS field (Flores 
et al. U999> . An attempt to detect a field population lensed by a 
massive cluster lens was also performed at 6.7 /mi (Altieri et al. 
1999 1, reaching a faint detection limit comparable to those of 
the HDF and the Lockman Hole. Most recently, a survey of the 
HDF-South (HDF-S; Oliver et al. 120021 has been reported with 
the same depth as that of the HDF. 

Here we present a very deep ISOCAM LW2 galaxy sur- 
vey in the SSA13 field. We first describe the observing strat- 
egy (Sect. |3 and the data processing (Sect.[3}. After evaluating 
the detection procedure by simulations (Sect. 0}, we present 
the catalog of detected sources (Sect. [3}. Finally we discuss 
the galaxy number counts and extragalactic background light 
(Sect.[6j. Other properties of the individual sources will be de- 
scribed elsewhere (Sato et al. 2003 1. 

2. Observation 

2.1. Field selection 

Small Selected Area 13 (SSA13), one of the Hawaii Deep 
Fields, was chosen as our deep 6.7/mi survey target. It is lo- 
cated at high galactic and ecliptic latitudes (b ~ 74° and 
P ~ 46°) and has a very low atomic hydrogen column den- 
sity of N(H I) = 1.5 X 10 20 cnr 2 (Cowie et al. U994t . Owing to 
these attractive properties for extragalactic studies, SSA13 has 
become one of the most-studied regions of the sky, comparably 
with the HDFs. It has photometric data obtained not only in 
the optical and near-infrared, but also at submillimeter and X- 
ray wavelengths (Cowie et al. [T996l Barger et al. [T9^S1 [T9^9l 
Mushotzky et al.|2000). A number of spectra were taken with 
the Low Resolution Imaging Spectrometer (LRIS) of the Keck 
telescope in a certain portion of the SSA13 field (Cowie et al. 
I1996 l. Because of the presence of a bright star (R = 10.4), the 
survey area was chosen to be the largest square region which 
is most distant from the star but inside the rectangular LRIS 
strip. The center of the region is (RA, Dec)= (13 h 12 m 26 s .O, 
42° 44' 24'.'8) in J2000. 

2.2. Observing strategy 

The survey observations were conducted with the mid-infrared 
camera ISOCAM (Cesarsky et al. H996b on board the Infrared 
Space Observatory ISO (Kessler et al. 1996 1. The sensitive 
broad band LW2 (5-8.5 fim) filter with a reference wavelength 
6.7//m (Moneti et al. 1 19971 1 was selected to detect radiation 
from stellar photospheres at high redshift. All the images were 
taken in the microscanning mode, CAM01, to minimize the 
flatfielding error (Siebenmorgen, Sauvage & Levine 1996 1. 
Total target-dedicated-time (TDT) for the observations was 23 
hours. 

Because of the visibility constraints for the SSA13 field, 
our observations were divided over six revolutions between 
revolution 560 and 614, 29 May 1997 and 22 Jul 1997 in UT 
(Table The observing time in each revolution was divided 
into two or three raster observations (Table |2). The observa- 



tions were executed essentially in two runs separated by a 1.5 
month interval, during which time the image plane of the space- 
craft rotated about 45 degrees (cf. roll angles in Table The 
main aim was to achieve better sampling of source profiles and 
also to maximize the overlapped area coverage. 

The surface density of bright sources suitable for image 
registration was expected to be low at 6.7 fim. Thus, it was al- 
most mandatory to fully utilize the unvignetted field-of-view 
for ISOCAM, 3 arcmin in diameter. We then selected the lens 
providing a 6 arcsec pixel field-of-view on the 32 x 32 pixel 
LW detector array. However, the use of 6 arcsec lens could in- 
troduce errors in the registration and hence degrade the depth of 
the map, because 6 arcsec pixels undersample the instrumental 
beam. At 6.7 /mi, the diffraction limited beam size is 2.3 arcsec 
FWHM (Full Width at Half Maximum) for the 60 cm ISO mir- 
ror. In order to suppress this undersampling effect, we applied 
several raster step sizes corresponding to non-integer multiples 
of the detector pixel size (Table|2ji. 

Because of the likely faintness of sources that could be used 
as references for the image registration, all the raster observa- 
tions in a single revolution were executed as a concatenated 
chain. If we had not adopted this concatenation setting, we 
could not align the CAM images and it would result in null re- 
sults. In the concatenated chain, the nominal automatic transi- 
tion to the standby mode between observations was suppressed. 
That transition, if permitted, would have introduced strong 
transients in the detector responsivity and random and system- 
atic displacements of the image on the detector (Siebenmorgen 
et al. 11996b . In fact, non-negligible displacements were actu- 
ally identified among revolution maps taken with different lens 
positions (TableQ. 

The integration time per exposure was set to T; nt = 20 sec. 
With shorter integration time, data would have been dominated 
by the readout noise of the detector because of very low back- 
ground emission at 6.7 /mi. The number of exposures per raster 
position was chosen to be N exp = 12, following a caveat for 
robust cosmic ray rejection (Siebenmorgen et al. 1996). The 
number of stabilization exposures was set to N sta b = 21 for the 
first raster in each revolution. The stabilization exposures were 
used as a buffer for the associated strong transients from the 
standby mode. Here the worst case was assumed, in which an 
observation would be preceded by a dark (zero flux) condition. 
For the subsequent rasters in any given revolution, the mini- 
mum stabilization value N sta b = 4 was used, since no change 
in background illumination would occur in the concatenated 
chains. The analog-to-digital converter gain was set to the max- 
imum x4 to cope with the low contrast of faint sources on the 
background. 

3. Image processing 

Mid-infrared data with ISOCAM suffer from cosmic rays (top 
panel of Fig.^. There is a huge number of spikes in the pixel 
histories and the signals slowly change even when the back- 
ground is constant. The latter is attributed to gradual changes 
or drifts in the detector responsivity due to the large number of 
cosmic ray impacts. Adopting the cosmic ray masking and the 
responsivity correction to deal with the above effects, we have 
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Table 1. Properties of the ISOCAM observations in the individual revolutions. Column 1 is the revolution number counted since 
ISO's launch. Column 2 is the total target-dedicated-time (TDT) in a single revolution. Column 3 is the roll angle of the satellite, 
counted counterclockwise from north to the raster N-direction. Column 4 is the solar aspect angle between the target and the 
Sun. Column 5 is the measured sky brightness derived with the sensitivity parameter of 2.32 ± 0.08 ADU/G/s/mJy for the LW2 
band (Blommaert 1998). ISOCAM signals in units of ADU/G/s (Analog-to-Digital Unit / Gain / second) can be converted to 
flux densities with the sensitivity parameter of the filter used. Column 6 is the lens position following the definition by Aussel 
et al. |r999l. Columns 7 and 8 are the measured offsets between the ISOCAM detector and the nominal center of the ISOCAM 
field-of-view determined by the telescope. They are shown in the raster M- and N-directions in units of pixels. These shifts were 
caused by the jitter in the lens wheel repositioning (Sect. l2~2l . 



Revolution 


Observing time 


Roll angle 


Solar aspect angle 


Sky 


Lens position 


Offset in M 


Offset in N 




(sec) 


(deg) 


(deg) 


(MJy sr 1 ) 




(pix) 


(pix) 


560 


11214 


315.2 


102.0 


2.90 


left 


0.95 


-0.24 


561 


11214 


314.4 


101.4 


2.94 


right 


-1.10 


-0.10 


610 


15315 


279.9 


69.2 


4.11 


left 


1.21 


0.16 


612 


15315 


278.4 


68.1 


4.21 


left 


0.90 


-0.12 


613 


15315 


277.6 


67.4 


4.28 


right 


-0.97 


0.19 


614 


15315 


276.7 


66.8 


4.19 


right 


-0.80 


0.14 



Table 2. Raster parameters used for the individual ISOCAM 
observations. Column 1 is the TDT number, the unique identi- 
fication number for ISO observations. Columns 2 and 3 are the 
numbers of raster points in the raster M- and N-directions. The 
raster directions were fixed in the spacecraft coordinate sys- 
tem. Columns 4 and 5 are the raster step sizes. Note the pixel 
field-of-view was set to 6 arcsec. 
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AN 








(arcsec) 


(arcsec) 


56000101 
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56000102 
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9 
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56100303 


6 


6 


7 
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56100304 


3 


2 


9 


9 


61000505 


6 


6 


7 


7 


61000506 


4 


4 


8 
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61000507 


3 


2 


9 
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61200108 


6 


6 


7 
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61200109 


4 


4 


8 


8 


61200110 


3 


2 


9 


9 


61300211 


6 


6 


7 


7 


61300212 


4 


4 


8 


8 


61300213 


3 


2 


9 


9 


61400314 


6 


6 


7 


7 


61400315 


4 


4 


8 


8 


61400316 


3 


2 


9 


9 



done the data processing following the scheme outlined below. 
Details of these steps are described in Appendix lAl They are 
coded as our own software. 

For each detector pixel, its raw signal raw can be written as 
the linear equation 



raw = a (/ s k y + /obj) + b, 



(D 



where a is the responsivity, b is the dark signal, and / S ky and 
/obj are fluxes from the sky and objects received by that pixel, 
respectively. The dark-subtracted signal dsub (top panel of 



Fig-D is then obtained as 

dsub — raw — b (2) 
= a (7 sky + 7 obj ). (3) 

The sky signal a / s k y can be extracted by masking the object sig- 
nal a 7 bj and also cosmic rays (Appendices IA. 3 .21 \AAl . This 
manipulation f(dsub) is essentially a linear interpolation of 
dsub signals applying the object and cosmic ray masks. The 
outcome reflects the detector response to the constant sky flux. 
We call it sresp (middle panel of Fig.Q, 

sresp = f(dsub) (4) 

= « ^sky, (5) 

and this will be used to correct the responsivity drifts. Then, 
the time-dependent responsivity a can be canceled out with the 
division of dsub by sresp, 

dsub I sresp = 1 + / bj// s ky (6) 
Thus, the fraction of the object flux ofrac can be given as 

ofrac = / ob j/7 s ky (7) 
= dsub I sresp - 1. (8) 

This technique for the responsivity correction works very well, 
resulting in a drastic reduction of noise (cf. bottom panel of 
Fig. ID- 

Although 7 s ky is constant in time, each detector pixel has 



a different / s k y due to the optical distortion (Aussel 1998) and 
vignetting in the LW array. Let the pixel receive a times larger 
sky flux than the fiducial pixel. Then, 



/s 



sky 



«4> 



(9) 



where /° ky is the sky flux for the fiducial pixel. The fiducial 
sky flux / ° kv is derived from the central part of the detector 
(Appendix lA.5l >. Then the ofrac data are transferred to the ob- 
ject flux ob j normalized to the fiducial pixel, 

obj = ofrac xJ°, (10) 
= /objM (11) 
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Fig. 1. An example of the responsivity correction to ISOCAM 
data taken with T; n t = 20 sec. The top panel shows the dark- 
subtracted data dsub = a(/ s ky + / bj)- The middle is sresp - 
a / s ky, the response to the blank sky after cosmic ray and ob- 
ject masking. The bottom shows ofrac - / bj//sky after cor- 
recting the responsivity drifts. For the sresp and ofrac panels, 
bad portions are masked out. The bad portions are cosmic rays 
themselves and portions where good linear interpolations are 
impossible with cosmic rays-masked and object-masked dsub 
data ( Appendices lA.3 . 31 [A.4> . 



Fig. 3. A search for the detection parameters (Sect. 14. 1> . 
The parameters for source detections are: the surface bright- 
ness threshold, the number of 0.6 arcsec sub-pixels above the 
threshold, and the FWHM of the convolution Gaussian filter. 
Diamonds denote a Gaussian FWHM of 3 sub-pixels, aster- 
isks 4 sub-pixels, squares 5 sub-pixels, and triangles 6 sub- 
pixels. A peak in N pos -N neg moves toward lower N pos with 
the larger FWHM, which corresponds to the suppression of 
noise detections. The lines indicate the detection reliability 
(N pos -N neg )/N pos = 100 % (solid line), 90, 80, 70, 60, 50, and 
40 % (dashed lines from the left to the right). 



Deriving the conversion factor for the fiducial pixel, we can 
obtain real object fluxes from the ob j data. 

After calculating the ob j data for all pixels, we made a map. 
To fully utilize the information of the raster positions and the 
optical distortion, we used sub-pixels of 0.6 arcsec to create a 
map from the obj data (Appendix lA.6> . After some iterations 
to refine the responsivity corrections with the cosmic ray and 
object masks and the offsets among revolutions, the final map 
was created. Fig.[2]shows the final map in the form of signal-to- 
noise ratio (S/N). The obj signals were co-added with weights 
of their inverse variances which were calculated by tracing the 
error propagation from the dsub stage (Appendix IA.2I >. Here 
we normalized standard deviations in sub-pixels in the coad- 
dition so that the central part of the histogram of the S/Ns 
in sub-pixels should resemble a Gaussian distribution with a 
sigma of one, because the noises in sub-pixels are correlated. 
There was no indication of a strong source confusion in the his- 
togram, supporting the use of a Gaussian fit to the histogram 
(Appendix lA.7> . 

4. Object extraction 

4.1. Determination of the detection parameters 

For extracting sources, the SExtractor program (Bertin & 
Arnouts 1996 1 was applied to Fig. [2 the final S/N map. The 
number of detected sources depends on three parameters: the 
surface brightness threshold, the number of 0.6 arcsec sub- 
pixels above the threshold, and the width of the Gaussian filter 
used to convolve the map. 



The detection parameters were determined by comparing 
numbers of detected sources on the S/N map with those on the 
negative S/N map (N pos and N neg , respectively). The negative 
S/N map was made by multiplying the S/N map by -1. Any 
detections made on the negative map are spurious, and no such 
sources should be detected (i.e. N neg = 0) where S/N is very 
high. The source extractions on the positive and negative maps 
were performed for more than 200 sets of the detection parame- 
ter values. The results are shown as a plot of N pos -N neg against 
Np 0S in Fig.0 

For a perfect Gaussian noise field with positive source pop- 
ulations, N pos -N neg should increase with N pos up to a certain 
level and then flatten out at some N pos , where the noise begins 
to be detected. In fact, Fig.[3]shows a declining trend where the 
flattening is expected. This indicates that the number of noise 
sources are somewhat larger in the negative side. One expla- 
nation may be due to our technique for responsivity correction 
(Sect. [3). A failure in masking an object or a cosmic ray will 
result in two dips in the neighboring data. Such dips are ac- 
tually identified around some bright sources in the final map 
(Fig. |2j. We believe that at least three dips are related to sources 
#54, #57, and #22 (Table[3}. Thus, their detection significance 
would be somewhat higher than the nominal values. 

We chose the detection parameters which gave the maxi- 
mum of N P os-N neg : 1.04 in the normalized S/N in sub-pixels 
for the surface brightness threshold, 119 sub-pixels above the 
threshold, and 5 sub-pixels for the width of the Gaussian. This 
parameter set resulted in N pos — 65 and N neg = 12. These val- 
ues give a nominal detection reliability of 82 % for the whole 
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1850_f2.jpg 

Fig. 2. The final map coadding all the data. North is top and east to the left in J2000. The center is (RA, Dec) = (13 h 12 m 26*0, 
42° 44' 2478). The map shows signal-to-noise ratio (S/N) per 0.6 arcsec sub-pixel. The signal is an average, weighted by the 
inverse of the assigned variance, and the noise is a normalized standard deviation (Sect. [3}. 



sample. However, the true detection reliability would be some- 
what higher as indicated by the decreasing trend in Fig. [3] 

4.2. Simulations of the detection procedure 

4.2.1. Detection completeness 

Monte Carlo simulations were performed to evaluate the com- 
pleteness of the source detection. For each simulation, a single 
test source was added to the real data, the S/N map (Fig. |2j. 
The profile of the test sources was fixed as a two dimensional 
Gaussian function. Its FWHM was determined by the real data; 
the intensity weighted average of the two smallest FWHMs 
among the ten most significant detections. It was 12 sub-pixels 
or 7.2 arcsec. 

Since we used the S/N map in the detection, we assigned 
S/Ns to the test sources with the random signal values. The 
noise was calculated for the 7.2 arcsec beam as a quadratic sum 
of the sub-pixel noises in that beam. This noise should be re- 
scaled, because the sub-pixel noises are not independent. The 
scaling factor was determined so that the histogram of the S/Ns 
in the 7.2 arcsec beam calculated over the entire map should 
be a Gaussian with a sigma of one. Thus, the noise includes all 
contributions from real sources (Barger et al.Qp99). The S/Ns 
were then determined with these signal and noise values. 

The simulations were done with almost 50,000 sets of ran- 
dom S/Ns and positions for test sources. As we used the real 
data, we set some conditions for detections; test sources should 
be detected within a distance of the FWHM and the excess of 
the measured fluxes over the assigned ones should be less than 
the assigned fluxes. In the case that a test source was added near 
a certain real source, some offset would be expected. However, 
if the offset becomes larger than the FWHM, this detection 
should be due to the pre-existing real source. Similarly, if the 
excess of the measured flux becomes larger than the assigned 
flux, such a detection should be related to pre-existing real 
source contributions. The result of the simulations is shown in 
Fig. |3 One sigma errors were estimated as Poissonian (Gehrels 
1986 1. The completeness of the detection is about 50% and 
80 % at assigned S/Ns of 2.5 and 4, respectively. 

4.2.2. Positional errors 

Fig-0presents the simulation results for the offsets of the mea- 
sured coordinates from the assigned ones. The offset is larger at 
lower S/N, and is larger than the expected values given by the 
formu la cr pos = FWHM/ V8ln2/(S/N) (dashed line; Condon 
119971 . A similar formula was used for a submillimeter map 
with a high source density of 1/11 sources per beam (Hughes 
et al. I1998i . Our map also has a high source density, but 1/30 
sources per beam. We suggest that the deviation from the theo- 




Input S/N 

Fig. 4. The completeness of the detection as a function of in- 
put S/N (Sect. l4.2.TV A single test source was added to the S/N 
map (Fig. |2j to evaluate the completeness of the detection pro- 
cedure. Almost 50,000 sets of S/Ns and positions were used. 
One sigma errors were estimated assuming Poisson statistics 
(Gehrels fT986l . 
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Fig. 5. The offsets of the output coordinates from the input 
ones. They are shown as a function of assigned S/Ns of the test 
sources. One sigma errors were estimated with eq. iA.H . Upper 
and lower sigma values were derived separately. The dashed 
line denotes the statistical formula for positional errors <x pos = 
FWHM/ V8ln2/(S/N) (CondonQjW). The size of sub-pixels 
in our map (0.6 arcsec) is also shown (dotted line). 



retical prediction is related to the high source density and also 
to the finite pixel size (dotted line) at high S/N. 
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4.2.3. Flux errors and biases 

Flux measurements always have errors. This effect can be well 
understood using the plot of assigned and measured fluxes nor- 
malized by their noise values in Fig.|6] In this plot of input and 
output S/Ns, we draw the unity line (dashed line) and one sigma 
deviation lines both for the positive and negative sides (dotted 
lines). The distribution of the simulation results on this plot is 
fairly well represented by these two dotted lines. Here we show 
simple moving averages and their associated one sigma errors 
calculated from eq. iA.H using solid circles with vertical er- 
ror bars. Large positive deviations from the unity line at low 
assigned S/Ns are mainly due to the incompleteness of detec- 
tions. Almost no sources are detected with measured fluxes less 
than two sigma. Even at relatively high S/Ns, there are slight 
positive deviations from the unity line. This could be due to 
effects of source confusion. 

It should be noted that this simple relation from input S/Ns 
to output S/Ns cannot be applied in the opposite direction: from 
measured fluxes to their true fluxes. In this case one must take 
into account the probability distribution of detected sources 
over this plot. Such a distribution can be derived by convolving 
the simulation results on this plot with the probability distri- 
bution along the horizontal axis (true sources), i.e., differential 
number counts, and also the completeness of detection (Fig.|4}. 

For the convolution calculation, we introduced fine 0.03 
dex meshes over the plot in accord with the steep slope of 
the number counts (Sect. 16 . II . All we needed here was rela- 
tive weights of the meshes along the horizontal axis at each 
vertical position. We assumed that the number counts were ex- 
pressed with a single power. In this case, the relative weights 
were determined by relative flux levels. Thus, we could use the 
S/N horizontal axis. We then applied the completeness values 
to these weights. Convolving the final weights with the simu- 
lation results, the transformation from measured S/Ns to statis- 
tically true S/Ns was derived and is shown as the red solid cir- 
cles with horizontal error bars in Fig. [6] Note that the assumed 
power of the differential number counts should be derived only 
after the correction of this flux bias. Thus, some iterations were 
needed here. 

Since the noises for the input and output S/Ns are measured 
on the noise map and are thus irrelevant to the simulations, this 
gives the relation between input and output fluxes. Because the 
slope in number counts is negative, fainter sources (less sig- 
nificant sources in this plot) are more numerous. As a result, 
the true flux bias (red circles) is larger than the simple relation 
(black circles) where the completeness is almost maintained. 
Beyond the point where the completeness begins to steeply de- 
cline to zero at the low input S/Ns or faint fluxes, the large 
deviations from the unity line seen in the simple relation begin 
to be suppressed. 

Flux biases for the detected sources are listed in Table [3J 
however, this correction is highly statistical. For example, if 
a particular source is known to be isolated, the correction de- 
rived here would be somewhat overestimated since the effect of 
source confusion would be smaller for that source. 




input S/N 

Fig. 6. The output S/N of the detected source as a function of 
input S/N. The individual simulation results are distributed al- 
most within the nominal one sigma limits (dotted lines) along 
the unity line (dashed line). Due to a combined effect of the 
incompleteness and the source confusion, the distribution is bi- 
ased toward higher output S/Ns. Means of the distribution and 
their one sigma errors are shown with solid circles with verti- 
cal error bars. For the correction of this flux biasing effect in 
the direction from output to input, the probability distribution 
over the diagram should be considered (Sect. l4~2~3l . It can be 
derived by convolving the simulation results with the differen- 
tial number counts and their completeness along the input S/N 
axis. The resulting relation from output S/Ns to input S/Ns is 
presented with the red solid circles with horizontal error bars. 

5. The catalog 

Table [3] lists the properties of all the 65 sources detected with 
the detection parameters described in Sect. 14.11 Here we ex- 
plain the absolute calibration of the coordinates and the fluxes. 

5.1. Astrometry - coordinates 

The coordinates in the catalog were calibrated with the USNO 
A2.0 catalog (Monet et al. I1998i extracted at the VizieR 
database (Ochsenbein 1997 1. Because very few LW2 sources 
were identified in the USNO catalog, the calibration was per- 
formed in two steps. The USNO catalog was first used to estab- 
lish the coordinate system in the / band image that was taken 
from the ground (Cowie et al. I1996i . The astrometric solution 
with 7 sources in the / band image shows a mean deviation of 
0725. The / band image was then compared with the LW2 im- 
age. Excluding confused sources, 12 sources in the LW2 band 
were used to derive the astrometric solution, giving a mean de- 
viation of 0775. Then the absolute astrometric error of the co- 
ordinate system in the LW2 image is 0779. The major cause of 
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Table 3. The catalog of ISOCAM LW2 sources in the SSA13 field. Column 1 is the identification number. Columns 2 and 3 are 
the Right Ascension in hours, minutes, seconds and the declination in degrees, arcmins, arcsecs in the J2000 system. Column 
4 is the total astrometric error (one sigma) in arcsec. It is a quadratic sum of the absolute astrometric error and the relative 
astrometric error calculated with 0.2 dex bins in Fig. [5] Column 5 is the normalized noise for the 7.2 arcsec diameter aperture 
at the source position. For its deviation, see Sect. 14.2.11 For the correction from instrumental values, a sensitivity parameter of 
2.32 ADU/G/s/mJy is used (Blommaert lT\138T l. Column 6 is the LW2 band (5-8.5/im) flux in juJy. This is the total flux with all 
the flux corrections applied: the flux bias (Sect. 14.2.31 . aperture, and processing corrections (Sect. 15 .21 . Taking account of the 
noise in all these corrections, their one sigma upper and lower errors are derived. SEDs are assumed to be flat (v/ v = constant). 
Column 7 is the flux bias defined as a ratio of true flux over measured one. Corresponding one sigma upper and lower limits are 
derived from Fig. [6] Column 8 is the completeness of the detection derived with 0.2 dex bins in Fig.|4] 
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this error is the combined effect of the significant undersam- 
pling of the instrumental beam with 6 arcsec pixels and the 
transients in the LW array (Sect. l5?2t . The total astrometric er- 
ror for an individual source is calculated as a quadratic sum of 
the absolute astrometric error and the relative astrometric error 
determined with Fig. [5] For sources with high completeness, 
this total astrometric error is conservative, because the effects 
of finite pixels have been accounted for twice. 



Fig-0shows the LW2 S/N contours superimposed on the / 
band image with the astrometric solution derived above. Many 
bright / band sources are easily identified with LW2 sources. 
We counted numbers of bright (/ < 23.7) sources for high 
and low completeness LW2 samples; there are 35 LW2 sources 
with completeness larger than 0.5 and 30 sources with lesser 
completeness. Within a distance of 1.5 times the total astro- 
metric error from the measured coordinates, 28 and 24 / band 
sources were found. The surface density of / < 23.7 sources 
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suggests that there would be 0.51 and 1.2 chance associations 
in the searched areas. Subtraction of these numbers gives rates 
of finding / < 23.7 sources as 79 and 76 % for the high and low 
completeness samples respectively. This small decrease in the 
rate indicates the high reliability of the catalog even for sources 
with low completeness. 



1850J7.jpg 

Fig. 7. LW2 S/N contours superimposed on the / band image. 
The image registration is performed using the astrometric solu- 
tion obtained in Sect. 15.11 The contour levels are set to 3, 4, 6, 
10 and 18 sigmas in sub-pixels (Fig.|2}. 



5.2. Photometry - flux 

The flux in the catalog is calculated statistically to be a total 
flux for a single point-like source. The total flux F tot ai is ob- 
tained with the formula F tota i = fbeam x 2.0 / 0.52. For each 



detected position, an aperture flux with 7.2 arcsec diameter is 
measured in units of ADU/G/s on the signal map. This flux is 
corrected using the flux bias determined in Sect. 14.2.31 This 
gives Fbeam- The two factors, 2.0 and 0.52, represent aperture 
and processing corrections which will be explained below. The 



Y. Sato et al.: A deep 6.7 fim survey in the SSA13 field with ISO 



9 



F tota i in units of ADU/G/s is then converted with a sensitivity 
parameter 2.32 ADU/G/s/mJy (Blommaert [l998t . 

The aperture correction was determined using the real data. 
Growth curves of the detected sources become noisy at diame- 
ters larger than 18 arcsec, but the dispersion was almost a min- 
imum at 7.2 arcsec. Thus, the aperture correction from a 7.2 
arcsec flux to an 18 arcsec flux was adopted. Similar to 7.2 
arcsec flux measurements (Sect. 14.2. 3l l. there would be flux bi- 
ases in 18 arcsec flux measurements. We therefore performed 
simulations for the 18 arcsec measurements similar to those for 
the 7.2 arcsec measurements, and corrected the flux bias in the 
18 arcsec measurements. The aperture correction factors were 
then obtained using the bias corrected 7.2 and 18 arcsec flux 
measurements for the 8 brightest sources, which resulted in a 
median of 2.0 with a one sigma error of 9.8 % (eq. iA.2\ ). At 
this stage, bias-corrected 18 arcsec fluxes were derived. 

Next, the processing correction factor was derived simulat- 
ing all the image processing steps by independent simulations 
different from those in Sect. 14.21 This processing correction is 
needed not only to translate the bias-corrected 18 arcsec fluxes 
to total ones but also to assess any effects in the image process- 
ing. Such processing effects are usually canceled out by a com- 
parison with observational data for external calibration sources 
that have been processed in the same manner. However, in this 
case no such data exist due to the specification of the ISOCAM 
observations. 

We therefore explain here the simulations used to derive 
the processing correction factor. At first, we needed to pre- 
pare simulated raw data. A total flux was assigned to a simu- 
lated source using the sensitivity parameter, and the source was 
put at a random position on the sky with the theoretical point 
spread function (PSF; Okumura 1998 ). It was then reprojected 
on the LW array at every exposure position using the satellite 
coordinates at every time step and the measured offsets among 
revolutions. The optical distortion coefficients of the LW ar- 
ray (Aussel [l99"8T l were also taken into account. The derived 
source signals were modified with the formula of the transient 
(eq. (6) in Coulais & Abergel l2000l using the measured back- 
ground flux (cf. Fig.|8j. Then, the measured responsivity drifts 
and dark signals based on the dark model were applied. Finally, 
we added the real ISOCAM raw data to these simulated raw 
data to simulate noise. This approach was mandatory because 
the noise properties in the ISOCAM data had not been properly 
understood. 

Next we applied the same image processing steps to these 
simulated raw data. The PSFs of simulated sources on the fi- 
nal map were found to be almost one arcsec sharper than those 
of real sources. This suggests that the technique for making 
simulated raw data was not perfect. One of the reasons for the 
sharper PSFs of the simulated sources could be errors in the 
measured offsets among revolutions. If the measured offsets 
were different from real values, PSFs of real sources would be- 
come broader artificially, while such effects cannot be involved 
in simulated sources. Actually, the magnitude of the differences 
in PSF widths is almost comparable to the total astrometric er- 
rors measured for real sources. There could be other possibil- 
ities such as errors in the optical distortion coefficients, theo- 




500 1000 1500 

Time (s) 



Fig. 8. An example of the transient effect of the ISOCAM LW 
array. The dashed line denotes the signal from an ideal detector 
free from the transient effect. A point source is assumed to pass 
over the center of the pixel with raster movement having a step 
size of 7 arcsec. The solid line represents the expected signal 
with the transient formula by Coulais & Abergel (2000). This 
figure assumes an instantaneous response of B = 0.51, a time 
constant of A = 560 ADU/G, a background flux of 8 ADU/G/s, 
and a source flux of 100/Jy. The vertical bar in the panel is a 
typical one sigma noise for each exposure. This implies that the 
transient correction for faint sources would only be performed 
statistically on the final product. 

retical PSFs which are too sharp, and ignoring cross-talk in the 
transient formula. 

In any case, the effects on the correction factor to the to- 
tal flux due to this discrepancy in the PSF shape will be min- 
imized for the large 18 arcsec aperture. We found that the 
bias-corrected 18 arcsec fluxes for the simulated sources cor- 
responded to 52+4 ' % of their total fluxes. There was no mean- 
ingful correlation between the correction factor and the input 
flux level. One of the main reasons for the large deviation from 
100 % may be our technique for responsivity correction. In this 
process, non-stabilized source signals are divided by responsiv- 
ities which are too high due to transient remnants (cf. Fig.|S|l. 

6. Number counts at 6.7//m 

6.1. Galaxy number counts 

To derive galaxy counts, we have to identify stars in the cata- 
log. We adopted SExtractor stellarity indices (Bertin & Arnouts 
1996 1 measured with the ground-based B and / band images 
(cf. Fig.0. The stellarity indices were calibrated with the HS T 
I band image having a partial coverage of the field (Cowie 
et al. U995a> . Two bright (> 100 (jJy) LW2 sources, #22 and 
#64 in the catalog, were then categorized as stars based on 
their high stellarity indices. At galactic latitude \b\ = 50°, 
Franceschini et al. ( I199U estimated that galaxies should dom- 
inate the 6.7 /im source counts below 90 fiJy. At this flux level, 
star counts should be 0.15 arcmin~ 2 , which is similar to our 
result of two stars in our survey area (16 arcmin 2 ). In the fol- 
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lowing discussion, all the sources in the catalog are assumed to 
be galaxies except #22 and #64. 

The number count is the number of sources in an unit solid 
angle with a certain flux constraint. It can be obtained from an 
observation as a sum of inverse of maximum solid angles in 
which the sources could be detected; 



(12) 



where N(S ) denotes a differential number count in a flux bin at 
flux S , i is a running number for detected sources in the corre- 
sponding flux bin, Q, is a maximum solid angle having a sensi- 
tivity to detect the source i. Integral number counts N(> S ) will 
be obtained by summing up N(S). In order to avoid unwanted 
effects due to negative sources, we restrict the summation to 
the significance level that could exclude negative sources. This 
level corresponds to a completeness of 56 %. Then the least sig- 
nificant source in the summation has a S/N of 4.3 in the map 
coaddition. 

Any incompleteness requires a calculation of effective solid 
angles. As the noise level is not uniform in our map, the effec- 
tive solid angle £2 should be obtained as an integral; 



Q(> S) 



I 



y(S /N)d£l(N), 



(13) 



where y denotes a completeness as a function of S /TV. The map 
is divided into small pieces d£l(N) giving the same noise level 
N. We convolve the non-uniform noise map with a 7.2 arcsec 
diameter aperture and relate them with the photometry of the 
detected sources using the aperture and processing correction 
factors (Sect. l5~2l . To go well with the restriction of the sum- 
mation, y is set to zero for y < 0.56. The derived effective solid 
angles £2(> S ) are shown in Fig.|5]as a function of flux. 

Then, integral galaxy number counts N(> S) at 6.7 /mi are 
obtained as listed in Table|4] The integral number counts climb 
up to 1.3 x 10 4 deg~ 2 at the faintest limit. The counts can be 
approximated with a single power -1.6 between 13/Jy and 
130/dy as 



N(>S) 



2.1!°* xlO 3 [deg- 2 ] 



40 [/dy] 



-1.6±0.3 



(14) 



As the integral counts are not independent of each other, we 
also list the differential galaxy number counts N(S) at 6.7 /mi 
in Tableg] 

The integral galaxy number counts are shown in Fig. EH 
with error bars. The vertical bars are derived with the Poisson 
statistics (Gehrels 1986). The horizontal bars were obtained 
applying mean logarithmic flux error factors for the detected 
sources counted in those flux bins. The error ranges in flux be- 
comes larger at fainter flux bins. 

6.2. Extragalactic background light 

Extragalactic background light (EBL) due to resolved galaxies 
can be derived as a first moment of the number counts: 



v7 v (5) = c/6.7/flnxJ]|^j, 




Flux S(6.7um) [uJy] 



(15) 



Fig. 9. Survey solid angles as a function of 6.7 /mi flux. 
Completeness of the detection is taken into account down to 
56 %. Below this limit, the completeness is set to zero. 



where vI v (S) is the 6.7 /mi EBL, Si is the 6.7 //m flux of the 
source i, and c is the speed of light. Integrating the extra- 
galactic contributions in the flux range of log S = 1.1-2.3 for 
our counts (Table |4j, we obtain the integrated 6.7 /jm light of 
0.56^|J [yhWirT 2 sr _1 . The upper and lower errors are derived 
as in Sect. 16.11 It should be noted that this integrated light is 
still a lower limit to the 6.7 /mi EBL, because the slope of the 
number counts is steeper than -1 even at the faint bins. So the 
6.7 //m EBL would become larger when integration to fainter 
flux levels becomes possible. 

Recently, TeV y-ray data were used to put constraints on 
the EBL in the infrared. This is based on electron-positron pair 
creation by photon collisions yxev + Tir — > e + + e~. Upper 
limit values on the infrared EBL are known to be sensitive to 
the assumed shape of the infrared EBL. If the infrared EBL 
had a local minimum around 5 /mi, the 5 cr upper limit could 
become as low as 2nWrrT 2 sr~ I (Renault et al. 2001 1. With 
a wide range of reasonable assumptions on the infrared EBL, 
a strict upper limit of 4.7nWirT 2 sr~' between 5 and 15 /mi 
was derived so far. From the galaxy modeling side, a value of 
2.4 nW iiT 2 sr~' at 6.7 /mi was predicted with the assumed red- 
shift of galaxy formation zf — 9 (Franceschini et al. 1991 1. 
Much more precise values of the EBL in this spectral range 
would put a strong constraint on the major epoch of star forma- 
tion in galaxies (Primack et al. U999> . 



6.3. Comparison with other 6.7 [im counts 

Our counts go somewhat deeper than previously reported 
counts at 6.7//m (Fig. 1 1 0i . In the overlapping flux range, 
our 6.7 /mi galaxy counts are almost consistent with the HDF 
counts (Oliver et al. 119971 Aussel et al. U999> . While the recent 
HDF-S counts (Oliver et al. 2002 not shown) are somewhat 
larger than the HDF and our counts, the field galaxy counts 
lensed by the massive cluster A2390 (Altieri et al. 19993 are 
about four times larger than ours. 



Y. Sato et al.: A deep 6.7 fim survey in the SSA13 field with ISO 



11 



Table 4. The galaxy number counts at 6.7 jum. N(> S ) are integral counts and N(S ) = AN(> S )/ A log S are differential counts. 
N pos (> S ) shows a number of the detected sources with a flux larger than S , excluding the two stars. N neg (> S ) is a number of the 
detected sources in the negative map. N pos (S) and N neg (S) are numbers of the detected sources in the corresponding flux bins. 
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The excess in the A2390 counts might be due to a failure in 
identification of cluster members. There may be errors in grav- 
itational lensing corrections in flux and solid angle. Multiple 
lensed images might be counted as individual field galaxies. 
These effects could become more significant since the target 
field is at the center of the cluster and the surveyed area is more 
than two times smaller than ours. Another reason might be con- 
tamination by stars. The SSA13 field and the HDFs are located 
at high galactic latitudes: b ~ 74° (SSA13), b ~ 55° (HDF), 
and b 49° (HDF-S). A2390 is located at low galactic lat- 
itude (b ~ -28°). At \b\ = 20°, stars are expected to domi- 
nate the counts even at a faint flux level of 3 fiJy (Franceschini 
et al. U991t . Altieri et al. ( 1999) stated that only two stars were 
identified among their 31 detections. The solid angles of the 
Altieri et al. (j 1999 1 counts might be underestimated by subtrac- 
tion of the pixel areas occupied with brighter sources. Aussel 
et al. dl999> did this and the nominal Aussel et al. ( 1999 1 count 
(green cross) is larger than the Oliver et al. ( 1997 1 counts in the 
HDF 

6.4. Comparison with models 

There are essentially two approaches to predict galaxy num- 
ber counts. One is a parametric evolution approach, which has 
been applied to IRAS 60 fim galaxy counts using IRAS 60 fim 
local luminosity function and its population mix. Normal spi- 
rals, starbursts, and AGNs are assumed to have non-evolving 
empirically determined SEDs but a parametric evolution in lu- 
minosity and/or density. After the observations with ISO and 
SCUBA, many authors have extended this modeling mainly to 
explain ISO 15 fim counts and SCUBA 850 /mi counts (Xu et 
al. |2t)03 and references therein). Some groups predict galaxy 
counts at 6.7 fim as well. Here we show the Rowan-Robinson 
(2001 1 model (dashed orange line in Fig. 1101. which could be 
marginally consistent with our counts. This model is known to 
represent well bright 6.7 fim counts. This will be due to the fol- 
lowing facts. In the local universe, the observing band of the 



ISOCAM LW2 (5-8.5/mi) filter matches the location of PAH 
and hot dust emission, which could be strong in star forming 
galaxies. Such low-redshift star-forming galaxies are known to 
be traced in the far-infrared, such as in the IRAS 60 fim band. 

Another approach is a stellar population synthesis method, 
which is used to model optical and near-infrared galaxy counts 
in searches for cosmological parameters. E/S0 and spiral 
galaxies are assumed and their local luminosity functions are 
adopted in the optical or near-infrared. Their SEDs evolve with 
time as their stellar populations do. Some include an AGN con- 
tribution as a parametric evolution population. Franceschini et 
al. ([1997 ) used this approach to model 6.7 fim counts (blue line 
in Fig. 1101 . Although the observational counts are associated 
with large errors, the consistency between this Franceschini 
et al. (1997 1 model and our results is remarkable. The steep 
slope comes from stellar emission of high redshift (z > 1) E/S0 
galaxies. At this redshift, the ISOCAM LW2 band will start 
to select near-infrared emission from stellar photospheres in 
galaxies, especially in E/S0 galaxies which are assumed to have 
formed most of their stars in the early epochs. Moreover, stel- 
lar SEDs short-ward of 6.7 /mi are steep up to 1 fim and their 
brightening at high redshifts is also expected. Because E/S0 
galaxies lack PAH and dust emission in the local universe, there 
is a negligible contribution from E/S0 galaxies in bright 6.7 fim 
counts. This affects additional contributions from E/S0 galaxies 
in this faint flux range. In contrast with its good agreement with 
the faint 6.7 fim counts, the Franceschini et al. ( 1997 ) model is 
known to underpredict bright 6.7 fim counts, probably due to 
insufficient evolution of dusty galaxies in the stellar population 
synthesis method. 

In this respect, a trial to combine these two approaches has 
been made by Xu et al. (2003 1. After constructing three para- 
metric evolution models, each of which is consistent with exist- 
ing dusty galaxy counts, they added E/S0 galaxy counts to each 
of them. The resulting counts have a slope consistent with our 
counts, but the counts themselves are much higher than ours, 
and in fact are more consistent with the A2390 counts. We 
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Fig. 10. The integral galaxy number counts at 6.7 fan. Our counts in the SSA13 field are marked with black filled circles with 
corresponding error bars (Sect. l6~Tl i. The initial results and the re-evaluation of the HDF data are shown as red hatches (Oliver 
et al.Q997 1 and green hatches and a cross (Aussel et al. H999i . The counts for field galaxies through the cluster gravitational lens 
A2390 are marked as cyan diamonds with their error bars (Altieri et al. 1999). Model predictions are overlaid for two models. 
The blue line represents the Franceschini et al. ( 1997) model and the dashed orange line is the Rowan-Robinson (2001 ) model. 



think that their population of starbursts should be merged with 
E/S0 galaxies at high redshifts. Xu et al. (2003) themselves re- 
port a synchronization of active evolutionary phases of E/SOs 
and starbursts at high redshifts (z > 1). A method of stellar 
population synthesis has now been extended to the far-infrared, 
and it shows that SEDs of E/S0 galaxies can resemble those of 
starbursts in their early evolutionary phases (Silva et al. I1998i . 
Interestingly, Pearson (2001 1 predicts excessive 6.1 /im counts 



similar to those of Xu et al. (2003) by adding a ULIG com- 
ponent to his parametric evolution model instead of an E/S0 
component. 

On the observational side, we find all submillimeter sources 
in the SSA13 field have SEDs consistent with high redshift 
(z > 1) galaxies having stellar masses comparable to those of 
typical E/SOs (Sato et al. 12002b . At the shorter wavelengths, a 
deficiency of high redshift (z > 1) elliptical galaxies was re- 
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ported from optical observations with spectroscopic and pho- 
tometric redshifts (Franceschini et al. I1998> . The number of 
morphologically identified elliptical galaxies at faint magni- 
tudes (/ > 24) was found to be smaller than the model predic- 
tions (Driver et al. 19981. Instead, the number of disturbed faint 
galaxies was reported to be larger and some of them were iden- 
tified with high redshift galaxies (Cowie et al. ll995al ll995bi. 
All these support the idea that elliptical galaxies undergo dusty 
starbursts at high redshifts. However, similar number count pre- 
dictions can be reproduced with a different set of SEDs and 
evolutionary parameters, especially in a parametric evolution 
approach as shown in the three parametric evolution models 
presented by Xu et al. (2003 1. Thus, a definite answer must 
await detailed investigations of individual faint 6.7 /jm galaxies 
detected in this survey (Sato et al. 2003 ). 



7. Summary 

A very deep mid-infrared survey has been conducted to investi- 
gate evolution of field galaxies at high redshift. The high galac- 
tic and high ecliptic latitude field SSA13 has been imaged with 
the ISOCAM's broad band filter at 6.7 pm. 

Utilizing quite high redundancy with a total observing time 
of 23 hours, we have developed a correction for responsivity 
drifts. We also introduced several types of masking and me- 
dian noise estimates for data affected by outliers (AppendixlAl. 
Many simulations have been carried out to properly account 
for the characteristics of the image processing and source de- 
tection. Within an area of 16 arcmin 2 , 65 mid-infrared sources 
down to 6 p Jy have been detected. For the central 7 arcmin 2 , an 
80 % completeness limit reaches 16/Jy. 

The galaxy number counts have been derived with statisti- 
cal corrections. The slope of the integral counts is - 1 .6 between 
13 //Jy and 130/Jy and the counts exceed 1 x 10 4 deg~ 2 in the 
faint limit. Our counts reach a limit about three times fainter 
than those in the HDFs. Although more detailed investigations 
are necessary, the steepness of the counts might indicate that 
stellar emission of evolving E/S0 galaxies has redshifted into 
the ISOCAM LW2 bands at this faint level. 
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Appendix A: Details of the image processing 

The image processing is summarized in Sect. [3] Here we de- 
scribe details of that processing; creation of raw, dsub, sresp, 
ofrac, and obj data, object masking, cosmic ray masking, and 
noise estimates with median absolute deviation. Map creation 
and noise investigation are also described. 

A.1. Dark subtraction 

Although all the processing steps in Sect.[5]are considered for 
one detector pixel, procedures in this subsection can be applied 
frame by frame. The raw signal frame raw in eq. Q was de- 
rived as a difference of ISOCAM image frames taken at the 
reset and the end of integration. Subtraction of the dark image 
frame b from the raw frame gave the dark-subtracted frame 
dsub. The dark image frames were created with the dark model 
invented and updated by Biviano et al. ( 1998) and Roman & 

Ott era. 

A.2. Median absolute deviation as a noise estimator 

The subsequent processing requires threshold values to identify 
objects and cosmic rays such as 3 sigma detection or 3 sigma 
rejection thresholds. In order to obtain such thresholds, one of 
the robust noise estimators - median absolute deviation - was 
adopted. Median absolute deviation t is defined here as 

t = med(abs(x; - med(jc,))) (A.l) 

for sample Xj. Detector noise, the sum of readout and photon 
noises, can be represented with the r, because the median op- 
erations in eq. jA.U are insensitive to outliers such as cosmic 
rays. The r was then scaled to give it the statistical meaning of 
the Gaussian sigma cr as 

cr T = t/t , (A.2) 

where to = 0.67448975. The scaling factor To is the t for a 
sample with a Gaussian distribution with cr = 1 . 

For the dsub data, <x T was calculated for each pixel in each 
revolution. With these <x T as starting values, all the error prop- 
agation was traced in the subsequent arithmetic operations. 

A.3. Responsivity correction 
A. 3.1 . Raster point masking 

When we are looking at blank sky, detector signals are domi- 
nated by the zodiacal emission. The zodiacal emission depends 
on the solar aspect angle (Table The temporal variation of 
the solar aspect angle is almost negligible during observations 
lasting less than several hours. Thus, the responsivity of the 
detector can be monitored by using blank sky as a reference. 

Because the responsivity drifts are major noise sources, cre- 
ating the response to sky sresp is the most essential part of the 
image processing. The sresp data were derived from the dark- 
subtracted data dsub. Because of the different sky levels among 
revolutions (Tabled, this operation was performed separately 
with the data taken in each individual revolution. 
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Fig. A.l. Schematic diagrams for the derivation of the response 
to the sky sresp. The top panel shows the derivation with raster 
point masking. The number of exposures per raster point is set 
to four (i.e. N exp = 4). Plus signs show dark-subtracted dsub 
signals of a detector pixel during five raster point observations. 
Vertical dotted lines show boundaries of raster points. A sresp 
value for exposure #5 (a red square on a red line) is derived with 
a line fitted to five dsub samples (red diamonds). In the bottom 
panel, object masks are applied. Four dsub samples at the sec- 
ond raster point are masked out (green asterisks). A sresp value 
for exposure #10 (a blue square on a blue line) is then derived 
correctly with a line fitted to nine unmasked dsub samples (blue 
diamonds). 



objects. The estimated noise is around 30 piy per pixel (1 <x) at 
this stage. Thus, almost no object could be identified in our ob j 
data (Sect.|52J. 



A low noise map must be derived in order to determine the 
location of the objects. All the obj data taken in a single revo- 
lution are then coadded ( Appendix IA.6i . The resulting map is 
called a revolution map. The noise found in the revolution map 
was more than ten times smaller due to coaddition of a few hun- 
dred ob j data in a revolution. At this stage, many bright sources 
could be identified with a certain threshold. Note that median 
values were used in the coaddition to suppress contamination 
by cosmic rays. 

The objects on the 2-dimensional (2-D) revolution map 
needed to be located in the 3-dimensional (3-D) data, because 
the sresp data are calculated with the 3-D dsub data. However, 
there is a complicated relation between sub-pixels in the 2-D 
revolution map and pixels in the 3-D data, due to raster slews, 
the optical distortion of ISOCAM and the use of sub-pixels in a 
revolution map (Appendix lA.6> . Thus, the projection of the 3-D 
ob j data to the 2-D revolution map was reversed. The resulting 
3-D data are called reversely-projected object robj data. Any 
pixel in the 3-D rob j data can be compared with that in the 3-D 
dsub data. This transformation from 2-D to 3-D is designed so 
that S/N in the 3-D rob j data is the same as in the 2-D revolu- 
tion map. A single S/N threshold was then used to mask pixels 
in the 3-D data looking at objects. 



The top panel of Fig. IA. II shows a schematic diagram for 
the sresp creation. A sresp value is to be assigned to each ex- 
posure. This value is evaluated via a linear interpolation of four 
or more dsub samples in adjacent raster points. For example, a 
sresp value for exposure #5 (a red square on a red line) is de- 
rived by fitting a straight line to three dsub samples at the first 
raster point (#1, 2, and 3) and two dsub samples at the third 
raster point (#8 and 9). Similarly, a sresp value for exposure 
#6 will be derived from five dsub samples (#2, 3, 8, 9, and 10). 
Because no dsub samples in its own raster point is used to de- 
rive a sresp value, we call this technique raster point masking. 

In the neighborhood of exposures with contributions from 
objects, however, the raster point masking gives higher sresp 
values than the correct ones. Such an example is shown for a 
sresp value for exposure #10 (a blue square on a blue line). In 
this setting, all the exposures at the second raster point should 
be masked out, because the detector is looking at a bright object 
at that raster point. With such masks, sresp values are then 
derived correctly. We will discuss such a procedure just after 
explaining object masks below. 

A.3.2. Object masks 

With the responsivity correction, division of the dsub data 
with the sresp data (eq. (jSJl), and the subsequent normalization 
(eq. ( I lOt ). the ob j data were obtained. They were calibrated to 
have zero flux for the sky. Thus, any positive features should be 



A. 3. 3. Raster point masking with object masks 

Now, we could apply the object masks to the dsub data. The 
raster point masking (Appendix IA.3.H was then repeated on 
the object-masked dsub data to refine the sresp data. 

The bottom panel of Fig. IA.1I shows a schematic dia- 
gram for the raster point masking with object masks. All the 
dsub data samples at the second raster point are object-masked 
(green asterisks). A sresp value for exposure #10 (a blue square 
on a blue line) is then derived from nine dsub samples #2, 3, 
12, 13, 14, 15, 16, 17, and 18 (blue diamonds). Note that a 
sresp value is estimated as the midpoint of the fitted line and 
at least two samples are required at one side of the fitted region 
to properly interpolate the responsivity drifts. 

In order to avoid inadequate interpolation, pixel history was 
regarded as discontinuous at positions having more than 3 cr 
deviation from the previous data. Near the discontinuities, the 
conditions of linear interpolation above can not be satisfied. 
Any exposures for which sresp values could not be defined 
were not used for the subsequent processing (Fig.[Q. 

The raster point masking with object masks was repeated 
three times. The object identification thresholds were de- 
creased; 5, 4, and 3 sigmas in this order. Note that no normal- 
ization of the noises in 6 arcsec pixels has been applied. At the 
final stage, the object masks occupy 3 % of the data. 
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Fig. A.2. A trial map constructed as if the data were taken with 
N eX p = I- The object masks were made from the final map 
(Fig. |2j with 2.5 sigma for the identification threshold param- 
eter. For the display, the noises in sub-pixels are normalized to 
correct the noise correlation in the map creation stage. 



A. 3. 4. Dwelling time and timescale of the responsivity 
drifts 

Some discussion is warranted concerning an appropriate value 
for the number of exposures per raster point; N exp = 12 for our 
case. Suppose one pixel looked at an object at the z™ raster 
point and it looked at sky at the adjacent raster points. The 
raster point masking routine calculates the sresp values for the 
i th raster point by interpolating the dsub data at the adjacent 
raster points. If the dsub data at the i th raster point are larger 
than the interpolated sresp data, this excess is recognized as an 
object. Because the responsivity is drifting, the noise caused in 
the interpolation becomes larger, if the time interval increases 
between the data at the I th raster point and the reference data at 
the adjacent raster points. 

To evaluate this effect, a trial map was created with an as- 
sumption that the data were taken as if N exp = 1 (Fig. lA.2l >. The 
object masks were created from the final map (Fig. [5) with an 
object identification threshold of 2.5 sigma. As expected, fluc- 
tuation in the background became roughly four times smaller 
than that in the final map. This indicates that smaller N exp (< 12) 
would be preferable to correct the responsivity drifts. Note that 
any objects below the identification threshold will never have 
appeared in this trial map. 



tamination from glitches, all the masks were canceled. New 
masks were then created from the ob j data. To preserve object 
fluxes in this step, a median of all exposures at each raster point 
was calculated as a reference value for masking. For example, 
a median of the ; th raster point was compared with all the data 
samples at this raster point. Any data samples which deviated 
from the threshold were masked out. 

The reference values for the first and second masking steps 
were derived from two dozen samples at most. Many points 
of sky within our field were observed more than 1000 times. 
Thus, in the third masking step, the low noise rob j data were 
utilized ( Appendix IA. 3. 2\ . All data in the obj data were com- 
pared with their respective rob j data. Any data were masked if 
they deviated from the threshold. 

Because we iterated the processing to refine object masks 
(Appendix lA.3.3> . the operation returned to the initial cosmic 
ray masking routine with the dsub data after the rob j masking 
step. In the second iteration and later, only the scale N = 2, 
instead of 3 and 11, was used to deglitch the dsub data. 

Different deglitching thresholds were set for cosmic ray im- 
pacts and their trailing remnants; typically 3 and 2.5 sigmas, 
respectively. These thresholds were adjusted to result in zero 
background in a map (eq. l|8}) and no loss in object fluxes. For 
discrimination of objects and cosmic rays, a standard devia- 
tion map and a median absolute deviation map (eq. (IA.2l i) were 
compared. Objects appear in both maps while cosmic rays ap- 
pear only in the standard deviation map. In the robj deglitch- 
ing, an additional margin for the thresholds was set to go well 
with the median nature of the robj data (Fruchter & Hook 
1998 1. At the final stage, cosmic ray masks occupy 23% of 
the data. 



A.4. Cosmic ray masks - deglitching 

The dsub data contain many spikes (top panel of Fig.[0. Such 
spikes, usually called glitches, were caused by cosmic ray im- 
pacts on the detector pixels. Glitches are short term phenom- 
ena, lasting a couple of exposures, and thus differ from the long 
term responsivity drifts. The deglitching algorithm was devel- 
oped to remove glitch effects in our data with T; nt = 20 sec. 
There are three masking steps based on dsub, obj, and robj 
data. 

Suppose that a glitch lasts N exposures in the dsub data. To 
remove this glitch with a time scale of N, the running median 
routine should be applied with a width of 2N + 1 . Any deviant 
data from the running medians were then regarded as glitches. 
The scale of glitches was typically N = 2; a cosmic ray impact 
and its trailing remnant. A new glitch sometimes occurred be- 
fore the preceding glitch had faded out. To take account of such 
overlaps, the scale was set to N = 3. Thus, 2N +1=7 was set 
for the width of the running median. This cosmic ray masking 
routine was repeated until no data deviated from the thresh- 
old. Then, an additional masking was applied with a scale of 
N = 11, equivalent to the number of exposures at a raster point 
(N exp = 12). 

This first masking step on the dsub data may mask bright 
objects. Thus, once the sresp data were derived without con- 



A.5. Sky brightness 

The fiducial sky flux T^L falling onto the fiducial pixels must 
be known to derive the object flux obj (eq. i!0\ ). The fiducial 
pixels are the central 12x12 detector pixels that have negligible 
effects from the optical distortion and the vignetting {a — 1 in 
eq. (|5)). The sky flux was dominated by the zodiacal light and 
the zodiacal light depends on the solar aspect angle (Table0. 
Thus, the sky flux must be measured for each revolution. The 
7^. was determined in two steps: (1) a median value of each 
pixel in the sresp data was calculated to minimize the effect of 
the responsivity drifts, and then (2) these median values of the 
central 12x12 pixels were averaged. The fiducial sky flux for 
each revolution is given in Tabled after the conversion to the 
sky brightness. 

A. 6. Map creation 

There exists an offset between the image on the detector and the 
actual satellite pointing. This is due to jitter in the lens wheel 
repositioning in the camera. The offsets due to the jitter vary 
from revolution to revolution. To evaluate these offsets, a rev- 
olution map was first created from the ob j data taken in each 
revolution. At first, we describe the creation of the revolution 
maps. 
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Information about the satellite position and the optical dis- 
tortion at the ISOCAM detector were used to coadd the obj 
data. To take account of the rotation and the distortion of the 
6 arcsec detector pixels on the sky, a map with 0.6 arcsec sub- 
pixels was created aligned with the J2000 system. The posi- 
tion and the orientation of the ISOCAM detector were obtained 
based on three angles defining the satellite status; RA, Dec, and 
Roll. Then, each detector pixel area on the sky map was de- 
fined with positions of the four corners of that pixel. The posi- 
tions were derived with the optical distortion correction depen- 
dent on the lens position (Aussel 1998 1. Finally, the distortion 
corrected pixel was shrunk to its 85 % (0.85 2 in area) and put 
into sub-pixels. This shrink operation recovers some informa- 
tion lost in the undersampling of the optical beam (Fruchter & 
Hook lT9m 

The pointing stability within T; nt = 20 sec was measured as 
0.34 arcsec (1 <x) on the average. Frames having displacements 
twice as large as this value were regarded as having been taken 
during slews of the satellite. Such frames were not used for the 
coaddition. The jitter of the lens wheel gave rise to a weakly 
illuminated area at some edges of the detector array. Pixels in 
such areas show 3 times larger noise in the error propagation 
calculations. They were masked out in the coaddition. The frac- 
tion of detector area affected by weak illumination depended on 
the lens position (Tabled; 3 % of a detector frame for the left 
position and 1 % for the right. 

After the revolution maps were obtained, positions of 
bright sources identified in different revolutions were com- 
pared. Based on the relative offsets of 19 sources, offsets among 
revolutions were derived (Table[0. With the correction of these 
offsets to the RA and Dec angles of the satellite, the final map 
(Fig. 13 was calculated in the same way to obtain revolution 
maps. Here all the ob j data in all the revolutions were used. 

Note that the relative offsets of the 19 sources had a mean 
deviation of 0.13 pixel, 0.8 arcsec. If such uncertainty is ac- 
ceptable, a priori jitter correction can be performed. Namely, 
the offset in the M-direction is 1 pixel for the left position and 
- 1 pixel for the right. The offset in the raster N-direction is 
(Table[D. 

A.7. Unfolding noise contributions 

The noise in the 7.2 arcsec beam was estimated to be typically 
1.1 /Jy/beam, with a use of a histogram for the correct abso- 
lute scaling of the noise. This noise includes detector, process- 
ing, and confusion noises. Here we tried to unfold these noise 
sources. 

For estimating a confusion noise, we adopted an evolution 
model by Franceschini et al. ( 1997 > because of the small de- 
viation of our counts from this model. For an extension to- 
ward fainter fluxes, we assumed a parameterized form for the 
Franceschini et al. model; 

2 N 



N(>S) = 



(*r + (*) G 



(A.3) 



with a set of parameters (ai,a2,So,No) = 
(0.70, 1.95, 20 fily, 6.4 x 10 3 deg~ 2 ). This formula repre- 
sents the model quite well at fluxes less than lOOyuJy. We 



generated test sources following this parameterized relation 
down to 0.05 /vJy and distributed them randomly over the 
sky. The fluctuation in the 7.2 arcsec beam over this map was 
estimated with a histogram, giving 0.49 yuJy/beam as a modal 
value. 

Using this simulated data set, we created raw data and all 
the image processing was executed as in the processing simu- 
lations in Sect. 15. 21 The resulting noise became 0.56 //Jy/beam. 
The increase in noise should come from a contribution of the 
processing noise, because we have not used real ISOCAM data 
in this case. By subtracting the confusion noise above quadrat - 
ically, the processing noise was estimated to be 0.28 /iJy/beam. 

As the total noise of the real ISOCAM data is 
1 . 1 //Jy/beam, the detector noise in the ISOCAM data should be 
0.91 yuJy/beam by a quadratical subtraction. It assures that our 
observations were dominated by the detector noise, i.e. readout 
and photon noises. 
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